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We present the hrst results from a new method for computing spacetimes representing corotating 
binary black holes in circular orbits. The method is based on the assumption of exact equilibrium. 
It uses the standard 3+1 decomposition of Einstein equations and conformal flatness approximation 
for the 3-metric. Contrary to previous numerical approaches to this problem, we do not solve only 
the constraint equations but rather a set of five equations for the lapse function, the conformal 
factor and the shift vector. The orbital velocity is unambiguously determined by imposing that, at 
infinity, the metric behaves like the Schwarzschild one, a requirement which is equivalent to the virial 
theorem. The numerical scheme has been implemented using multi-domain spectral methods and 
I passed numerous tests. A sequence of corotating black holes of equal mass is calculated. Defining 

the sequence by requiring that the ADM mass decrease is equal to the angular momentum decrease 
multiplied by the orbital angular velocity, it is found that the area of the apparent horizons is 
^ , constant along the sequence. We also find a turning point in the ADM mass and angular momentum 

_ O ■ curves, which may be interpreted as an innermost stable circular orbit (ISCO). The values of the 

' global quantities at the ISCO, especially the orbital velocity, are in much better agreement with 

those from third post-Newtonian calculations than with those resulting from previous numerical 
approaches. 
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I. INTRODUCTION 

O . Motivated by the construction of several gravitational wave detectors (LIGO, GEO600, TAMA300 and VIRGO) 
great efforts have been conducted in the past years to compute the waves generated by binary black holes. We 
r presented in Ref. (hereafter Paper I) a new method for getting quasi-stationary spacetimes representing binary 
black holes in circular orbits. See also Paper I for a review on issues and previous works in this field. 
The basic approximation is to assume the existence of an helical Killing vector : 



dto aipo 

where d/dto (resp. d/dipo) is a timelike (resp. spacelike) vector which coincides asymptotically with the time 
coordinate (resp. azimuthal coordinate) vector of an asymptotically inertial observer. Basically, it means that the 
two black holes are on circular orbits with orbital velocity fl This is of course not exact because the emission of 
gravitational waves will cause the two holes to spiral toward each other. But this is a valid approximation as long as 
the time-scale of the gravitational radiation is much longer than the orbital period, which should be true, at least for 
large separations. The existence of £ enables us to get rid of any time evolution. 

We use the standard 3+1 decomposition of the Einstein equations We restrict ourselves to a space metric that 
is conformally flat, i.e. of the form : 

7 = *'f, (2) 
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where is a scalar field and f denotes the flat S-metric Let us mention that the exact spacetime should diff'er 
from conformal flatness and that this assumption is only introduced for simplification and should be removed from 
later works. However it is important to note that it is consistent with the existence of the helical Killing vector and 
the assumption of asymptotic flatness. The ten Einstein equations then reduce to five equations, one for the lapse 
function N, one for the conformal factor 4" and three for the shift vector (3 (see Paper I for derivation) : 



A7V = N^^AijA'^ - 2Dj In 'i'D^N 
A/?* + \&Dj(}i ^ 2A'^ {DjN ~ 6NDj In*) 



Aij A 



(3) 
(4) 

(5) 



where Di denotes covariant derivative associated with f and A :— D^D^ the ordinary Laplace operator. A^^ is the 
reduced extrinsic curvature tensor related to K^^ by A^^ := '^'^K^^ and given by 

A^'=^{LI3f, (6) 

{Lpy-' denoting the conformal Killing operator applied to the shift vector 

iLf3y' &f3' + - ^Dk^T^. (7) 

Equations (||), (||) and (||) are a set of five strongly elliptic equations that are coupled. To solve such a system, we 
must impose boundary conditions. To recover the Minkowski spacetime at spatial infinity, i.e. asymptotical flatness, 
the flelds must have the following behaviors : 



N ■ 



1 when 



fl— — when r 

difo 
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oo 



4' 



1 when 



(8) 
(9) 
(10) 



As we wish to obtain solutions representing two black holes and not Minkowski spacetime, we must impose a non- 
trivial spacetime topology. In Paper I, we define the topology to be that of the real line M times the 3-dimensional 
Misner-Lindquist manifold |^,|| ; this defines two throats, being two disjointed spheres 5*1 and S2 of radii ai and 02, 
centered on points (xi, 0, 0) and (0:2, 0, 0) (such that \xi — X2\ > 01+02). Following Misner Lindquist Kulkarny 
et al. §, Cook et al. |-0 and others we demand that the two sheets of the Misner-Lindquist manifold are 

isometric. Moreover we choose the lapse function N to be antisymmetric with respect to this isometry. We solve the 
Einstein equations only for the "upper" sheet, i.e. only for the space exterior to the throats, with boundary conditions 
given by : 



= 



and 
and 



N\ 



S2 







/3 







dri 2ri 



and 



S2 

dr2 2r2 



0, 



(11) 
(12) 

(13) 



S2 



where ri and r2 are the radial coordinates associated with spheres Si and 5*2. Equations (0) reflect the antisymmetry 
of the lapse function TV. The boundary conditions for the shift vector, given by Eqs. (|12|), represent two black holes 
in corotation (rotation synchronized with the orbital motion), which is the only case studied in this paper. Those 
boundary conditions should be easily changed to represent other states of rotation (like irrotation) . Equations ( |l^ ) 
come from the isometry solely. 

The orbital velocity 17 only appears in the boundary condition for the shift (see Eq. (H)). Equations (||), (||) and 
(^ can be solved for any value of SI. So we need an extra condition to fix the right value for Vl. This is done by 
imposing that, at spatial infinity, the metric behaves like a Schwarzschild metric, i.e. by imposing that has no 

monopolar term in l/r : 



^^N - 1 



— when r 



(14) 
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In other words, is chosen so that the ADM and the "Komar-Hke" masses coincides, those masses being given by 

Madm - -7^ / D'-^dS, (15) 

Joe 

MKomar = ^ / D'NdS,. (16) 



As shown in |13| and in Paper I this is closely linked to the virial theorem for stationary spacetimes. We will see 
later that this uniquely determines the orbital velocity, and that this velocity tends to the Keplerian one at large 
separation. 

This paper is organized as follows. Sec. p is dedicated to the presentation of the numerical scheme, that is based on 
multi-domain spectral methods. In Sec. |ll we present some tests passed by the code, which encompass comparison 



with the Schwarzschild and Kerr black hole and the Misner-Lindquist solution \^^- In Sec. IV we present results 
about a sequence of binary black holes in circular orbits. In particular we locate the innermost stable circular orbit and 
compare its location with other works. Sec. ^ is concerned with extension of this work, for getting more complicated 
and more realistic results. 

II. NUMERICAL TREATMENTS 
A. Multi-domain spectral methods 

The numerical treatments used to solve the elliptic equations presented above is based on the same methods that 
we already successfully applied to binary neutron stars |14| . The sources of the equations being mainly concentrated 
around each hole we use two sets of polar coordinates centered around each throat (see Sec. ^. Note however that 
the tensorial basis of decomposition is a Cartesian one. For example, a vector field V will be given by its components 
on the Cartesian basis {Vx, Vy,Vz) but each component is a function of the polar coordinates (r, 6, ip) with respect to 
the center of one hole or the other. 

We use spectral methods to solve the elliptic equations presented in Sec. | ; the fields are given by their expansion 
onto some basis functions. Mainly, we use expansion on spherical harmonics with respect to the angles (6, ip) and 
Chebyshev polynomials for the radial coordinate. Let us mention that there exists two equivalent descriptions : a 
function can be given in the coefficient space, i.e. by the coefficients of its spectral expansion, or in the configuration 
space by specifying its value at some collocation points |l^ . 

The sources of the elliptic equations being non-compactly supported, we must use a computational domain extending 
to infinity. This is done by dividing space into several types of domains : 

• a kernel, a sphere containing the origin of the polar coordinates centered on one of the throats. 

• several spherical shells extending to finite radius. 

• a compactified domain extending to infinity by the use of the computational coordinate u — K 

This technique enables us to choose the basis function so that the fields are regular everywhere, especially on the 
rotation axis and to impose exact boundary conditions at infinity. This has been presented with more details elsewhere 
H^6|- [l8|. Note that since the last domain extends to spacelike infinity, the surface integrals defining global quantities. 



such as (15) and (|T6|), can be computed without any approximation. This contrasts with other numerical methods 
based on finite domains (cf. e.g. Ref. |J9| and Fig. 1 therein). As two different sets of coordinates are used, one 
centered on each hole, we are left with two computational domains of this type, each describing all space and so 
overlapping. 

The sources of the equations being concentrated around the two throats, we wish to split the total equations (||), (^) 
and (H) into two parts, each being centered mainly around each hole and solved using the associated polar coordinates 
set. So an equation of the type AF — G will be split into 

A^^i = Gi (17) 
AF2 = G2, (18) 

with F = F1+F2 and G = Gi +G2. Ga is constructed to be mainly concentrated around hole a, and so well described 
by polar coordinates around this hole. Therefore, the solved equations are : 
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AiV, = N-^^A.jA"^ - ^D'-^aDjN (19) 



A(3l + \&D,(3i = 2A^^ {D,Na - 6^ (20) 

A*a--^i.,ijf, (21) 

where the values with no index represent the total values and the values with index a represent the values "mostly" 
generated by hole a (a = 1 or 2). For example, we have DiN — DiNi + DiN2, DiNa being concentrated around hole 
a. Doing so, the physical equations and sources are given by the sum of equations (|2^), ( |20|) and ( |l9|) for a = 1 and 
a = 2. For more details about such a splitting of the equations into two parts we refer to |14[]. 

B. Elliptic equations solvers 

1. Scalar Poisson equation solver with boundary condition on a single throat 

Using spectral methods with spherical harmonics, the resolution of the scalar Poisson equation reduces to the 
inversion of banded matrices. We already presented in details in ||l^,|l^ the methods to solve such equations in all 
space, imposing regularity at the origin and exact boundary condition at infinity. In the case of black holes we wish 
to replace the regularity at the origin by boundary conditions on the spheres Si and 5*2 and to solve only for the part 
of space exterior to those spheres. In Ref. ||l^ we have shown that, for each couple of indices {I, m) of a particular 
spherical harmonic, we can calculate one particular solution in each domain, two homogeneous solutions in the shells 
and only one in the kernel (due to regularity) and one in the external domain (due to boundary condition at infinity) . 
The next step was to determine the coefficients of the homogeneous solutions by imposing that the global solution is 
at the boundaries between the different domains. 

In the case of a single throat S, the boundary condition is given by a function of the angles solely, i.e. B {9, ip). One 
wishes to impose that the solution or its radial derivative is equal to B on the sphere which corresponds respectively 
to a Dirichlet or a Neumann problem. We choose the kernel so that its spherical boundary coincides with the throat. 
So we do not solve in the kernel with represents the interior of the sphere. B is expanded in spherical harmonics and 
for each couple {l,rn), we use one of the homogeneous solution in the first shell to fulfill the Dirichlet or Neumann 
boundary condition. After that we are left with one particular solution in every domain, one homogeneous solution 
in the innermost shell and in the external domain and two in the other shells. The situation is exactly the same as 
when a solution was sought in all space and the coefficients of the remaining homogeneous solutions are chosen to 



maintain continuity of the solution and of its first derivative. So the generalization of the scheme presented in 18 
is straightforward and enables us to solve either the Dirichlet or Neumann problem, with any boundary condition 
imposed on the throat. 

2. Vectorial Poisson equation solver with boundary condition on a single throat 

We presented extensively two different schemes to solve the vectorial Poisson equation (^ in all space in ||l^ (the 
Oohara-Nakamura j2^] and Shibata |^ schemes). We present here an extension of the so-called Oohara-Nakamura 
scheme to impose boundary condition a throat and to solve only for the exterior part of space. The Shibata scheme 
has not been chosen because, the solution being constructed from auxiliary quantities, it is not obvious at all to 
impose boundary conditions on it. This is not the case with the Oohara-Nakamura scheme where the final solution is 
calculated directly as the solution of three scalar Poisson equations. More precisely the solution of (cf. Eq. ^0|) 

A/3* + XD'Djl3^ = V' (A 7^ -1) (22) 

is found by solving the set of three scalar Poisson equations 

A(3' ^V' ~ XD'x, (23) 

where x is solution of 

Ax=j^D,V\ (24) 
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Let us mention that this scheme should only be used with a source V that is continuous. We use the scalar Poisson 
equation solvers with boundary condition previously described to solve for each Cartesian component of (|23|) with the 



appropriate boundary conditions. But let us recall (see 1 18 ) that the Oohara-Nakamura scheme is only applicable if 



X = (25) 

and that it only ensures that 

A (x - A/30 = 0. (26) 
One can easily show that ( p6| ) implies (^5|) if and only if 

x\s=D^(3'\s' (27) 



which is the boundary condition we must impose during the resolution of (24) to use this scheme. Let us mention 
that X being calculated before (3, we must use some iterative procedure. We first solve (^4|) with an initial guess of 
the boundary condition and then determine /? by solving (p3|). Using that value, we can determine a new boundary 
condition for x, using (^^, and so a new (3. This procedure is repeated until it has sufficiently converged. The 
obtained (3 is then solution of the vectorial Poisson equation with either a Dirichlet or Neumann type boundary 
condition on the sphere S. 



3. Elliptic solvers with boundary conditions on two throats 



In order to illustrate how boundary conditions arc put on the two spheres Si and ^2, let us concentrate on the 
Dirichlet problem for the scalar Poisson equation. One wishes to solve 



AF = G, 



with the boundary conditions 



^1, 
F\ 



S2 



Bi (01, (^i) 

B2 (02,^2) , 



where Bi and B2 are arbitrary functions. As explained in Sec. II A, the total equation is split into two parts 



AFi = Gi 
AF2 = G2, 



(28) 



(29) 
(30) 



(31) 
(32) 



the equation labeled a = 1 or 2, being solved on the grid centered around hole a so that the sphere Sa coincides with 
the innermost boundary of the first shell. 

During the first step we solve Eqs. (31) and (p3) with the boundary conditions 



(33) 
(34) 



by means of the scalar Poisson equation solver described in Sec. 
does not fulfill the boundary conditions (p9|)-(p 



[IBl 



boundary condition (34) by = B2 — Fi 



Doing so, the total solution F — Fi + F2 
|30| ). So we calculate the values of Fi on the sphere S2 and modify the 
The same modification is done with the boundary condition (^3|) . Then 
we solve once again for Fi and F2. The whole procedure is repeated until a sufficient convergence is achieved. So 
we are left with a function F which is solution of the Poisson equation (|8|) and which fulfills a given Dirichlet-type 
boundary condition on two spheres (p^)-(pO|). 

The same thing can be done for the Neumann problem by modifying the boundary conditions using the radial 
derivatives of the functions Fa. The same technique is applied for the vectorial Poisson equation. Let us mention that 
the iteration on the boundary conditions for /?, resulting from the presence of the two spheres, is done at the same 
time than the one on the quantity x resulting from to the Oohara-Nakamura scheme (see Sec. [IB 2). 
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4- Filling the interior of the throats 



As seen in the previous section, we can solve elliptic equations with various boundary conditions in all the space 
exterior to two non- intersecting spheres Si and 5*2. But a problem arises from the iterative nature of the total 
numerical procedure. Suppose that after a particular step the lapse N = Ni + N2 has been calculated by means of 
the two Poisson equations ([T9|). From the very procedure of the elliptic solvers, iVi (resp. N2) is known everywhere 
outside sphere 5*1 (resp. 5*2). If the next equation to be solved is the one for the shift vector split like (^o|), N appears 
in the source term. We need to know the source everywhere outside the associated sphere Sa (a = 1, 2) which includes 
the interior of the other sphere. So we must construct fields that are known in the all space. After each resolution, 
the fields are filled as smoothly as possible inside the associated sphere. In our example, after the resolution of (|is|), 
A^i and N2 are filled inside the spheres, so that the total function N is known everywhere. 

The filling is performed, for each spherical harmonic {I, m), by the following radial function : 

• (Sr"* - 2rS) (a + /Jr^) if I is even, 

• (Sr-* - 2r^) {ar + /3r^) if I is odd, 

where the coefficients a and /? are calculated so that the function is across the sphere Sa- The multiplication by 
the polynomial (Sr"* — 2r^) ensures that the function is rather regular at the origin. Of course this choice of filling is 
not unique and the final result should be independent of the filling procedure, the fields outside the spheres depending 
only on the boundary conditions on those spheres. The choice of filling may only change the convergence of the 
numerical scheme. Let us stress that even if the fields are known, regular and everywhere, they have a physical 
meaning only outside the throats. The filling is only introduced for numerical purposes. 



C. Treatment of the extrinsic curvature tensor 



1. Regularization of the shift 



When one imposes corotation for the two black holes, that is a vanishing shift vector on the throats, isometry 
conditions (59), (60) and (61) of Paper I are trivially fulfilled. Unfortunately this is not the case for (62) and (63) of 
Paper I. We must find a way to impose that 



dr 



= and 



Si 



dr 



= 



(35) 



Si 



in order to get a truly isometric solution. 

Another problem comes from the computation of the reduced extrinsic curvature tensor A*-' by means of (|^). 
Because of division by iV = on Si and S2, we must impose that 



Si 



and 



S2 



0, 



(36) 



so that the extrinsic curvature tensor is regular everywhere. Because of the rigidity conditions ( |12D and for a truly 
isometric solution verifying (|35|), the regularity conditions (|3^) are fulfilled if and only if 



a/3'- 



dri 



and 



a/3'- 



Si 



dro 



0. 



(37) 



So, to get a truly isometric and regular solution, both the value and the radial derivative of /3 must be zero on the 
throats : 



and 



a^ 

dr 



0. 



(38) 



s. 



But when solving Eq. (|j), one can only impose the value at infinity and one of those two conditions, i.e. we can 
only solve for the Dirichlct or Neumann problem, not for both. We choose to solve the equation (^ for the Dirichlet 
boundary condition : /3 = on both spheres. Doing so, the regularity conditions (pTh, as well as the remaining 
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isometry conditions 
enforce ( |37| ) and (|35 



|3^), are not necessarily fulfilled. After each step we must modify the obtained shift vector to 
The part of the shift generated by the hole 1 is modified by 



r-'c' 



{R - r,f (n - ai) d 



old 



dri 



(39) 
(40) 



where ri is the radial coordinate associated with hole 1, ai the radius of the throat and R an arbitrary radius, typically 
R — 2ai. The correction procedure is only applied for ai < ri < R. Let us mention that the function of ri in front 
of d /dri in Eq. (^) has been chosen so that it maintains the value of the shift vector on the sphere 1 and its 

continuity (C^ fimction). The same operation is done for the other hole. After regularization, the shift vector satisfies 
(i) the rigidity condition (p^, (ii) the isometry conditions (|35|), and (iii) the condition ( p7| ) ensuring the regularity of 
the extrinsic curvature, but it violates slightly the momentum constraint (^). 
As seen in Paper I, the regularity is a consequence of the equation 



A/3' 



-6/3' A In*. 



(41) 



Because this equation is not part of the system we choose to solve, w e do no t expect that the correction function 
is exactly zero at the end of a computation. But we will verify in Sec. IVB 1 that it is only a small fraction of the 
shift vector (less than 10^"^), fraction which represents the deviation from Eq. (^ ) (sec also Ref. for an extended 
discussion). Moreover, we will see in Sec. [II B that /3cor converges to zero for a single rotating black hole. 



2. Computation of the extrinsic curvature tensor 

Using the regularized shift vector presented above, we can compute the tensor (i/3)'"' , which is zero on both throats. 
To calculate the tensor A^^ one must divide it by the lapse function which also vanishes on both throats. Near the 
throat 1 , has the following behavior 



^In^ai = (ri -ai)ni , 



(42) 



where ni is non zero on throat 1 (this supposes that ri = ai is only a single pole of A^, which turns out to be true, 
dN/dri representing the "surface gravity" of black hole 1). We can compute rti, using an operator that acts in the 
coefficient space of A^ and divides it by (ri — ai). The same operation is done with 



(ri - ai) ll^ 



(43) 



The divisions are also done on the second throat. To compute the extrinsic curvature tensor in all space we use 

• A'^ — V-^ / (2ni) in the first shell around throat 1 

• A^^ — I2 I (2n2) in the first shell around throat 2 

• = [Lpf^ I {2N) in aU other regions. 

This procedure enables us to compute the extrinsic curvature tensor everywhere, without any problem that could 
arise from a division by zero. 



3. Splitting of the extrinsic curvature tensor 

In the split equations (|l^) and (|2ll), the term A^ appears. This term represents the part of the total extrinsic 
curvature tensor generated mostly by hole 1 so that the total tensor is given by 

A'^^A{^+Al^. (44) 

For the binary neutron stars treated in those split quantities were constructed by setting Af = {LPif'' / (2N). 
Such a construction is not applicable in the case of black holes. Indeed, only the total shift vector is such that 
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{LPY'' = on the throats and not the spht shifts /3i and [32- If such a construction were apphed the quantity A^^ 
would be divergent due to division by A'^ = on the throats. The computation presented in the previous section gets 
rid of such divergences but enables us to calculate only the total A^^ . 
The construction of A^^ and A^^ is then obtained by 

A^f = A'^Hi (45) 
i^^ = A''H2, (46) 

where Hi and H2 are smooth functions such that Hi + H2 — 1 everywhere. We also want Hi (resp. H2) to be close 
to one near hole 1 (resp. 2) and close to zero near hole 2 (resp. 1), so that A^ (resp. A2 ) is mostly concentrated 
around hole 1 (resp. 2). So, we define Hi by 

• 1 if ri < Riat 

• \ [C0S2 (^/2 (n - iJint) / (i?cxt - iiint)) + l] if i?int < < i?cxt 

• if r2 < i?int 

• \ Sin^ (7r/2 (ri - i?int) / (i?oxt - i?int)) if ^int <T2< i?cxt 

• i if ri > i?oxt and r2 > i?cxt, 

where ri (resp. r2) is the radial coordinate associated with throat 1 (resp. 2). The radii i?int and i?cxt are computa- 
tional parameters, chosen so that the different cases presented above are exclusive. Typically, we choose -Rint = d/6 
and i?cxt = d/2, where d is the coordinate distance between the centers of the throats. H2 is obtained by permutation 
of indices 1 and 2. 



D. Numerical implementation 

The numerical code implementing the method described above is written in LORENE (Langage Objet pour la 
RElativite NumeriquE), which is a CH — h based language for numerical relativity developed by our group. A typical 
run uses 12 domains (6 centered on each black hole) and NrXNgx = 33 x 21 x 20 (resp. NrXNgX = 21 x 17 x 16) 
coefficients in each domain in high resolution (resp. low resolution). For each value of fl, a typical calculation takes 50 
steps. To determine the right value of the angular velocity, by means of a secant method, it takes usually 5 different 
calculations with different values of fl. The associated time to calculate one configuration is approximatively 72 hours 
(resp. 36 hours) for the high resolution (resp. for the low resolution) on one CPU of a SGI Origin200 computer 
(MIPS RIOOOO processor at 180 MHz). The corresponding memory requirement is 700 MB (resp. 300 MB) for the 
high resolution (resp. low resolution). 



III. TESTS PASSED BY THE NUMERICAL SCHEMES 



A. Schwarzschild black hole 



In this section we solve Eqs. (||) and (||), with boundary conditions (|T^) and (|T^) on a single throat S. The 
behaviors at infinity are given by Eqs. (||) and (^^. In this particular case, the shift vector P is set to zero, so 
that A^^ vanishes. This represents a single, static black hole, and we expect to recover the Schwarzschild solution in 
isotropic coordinates. 

The computation has been conducted with a initial guess far from the expected result. More precisely, we began 
the computation by setting iV = 1 and = 1 everywhere. Equations (||) and (||) are then solved by iteration. Let us 
mention that the boundary condition on the conformal factor, given by (^j), is obtained by iteration. At each step 
we impose 



9?- 



,7-1 



2r 



(47) 



where is the conformal factor at the current step and 5''^ ^ at the previous one. 
Before beginning a new step, some relaxation is performed on the fields by 
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Q'' ^ \Q' + (1 - \)Q 



j-i 



(48) 



where < A < 1 is the relaxation parameter, typicahy A = 0.5. Q stands for any of the fields for which we solve an 
equation {N and ^' solely for the static case). 

The iteration is stopped when the relative difference between the lapse obtained at two consecutive steps is smaller 
than the threshold SN — 10~^'^. The computation has been performed with various number of collocation points and 
with two shells. All the errors are estimated by the infinite norm of the difference. 
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FIG. 1. Relative difference between the calculated and the analytical lapse A'^ with respect to the number of radial spectral 
coefficients for the Schwarzschild black hole. The circles denote the error in the innermost shell, the squares that in the other 
shell and the diamonds that in the external domain. 
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FIG. 2. Same as Fig. |l| for the conformal factor 'P. 

Figures |] and || show a extremely good agreement with the exact analytical solution. The saturation level of 
approximatively 10~^^ is due to the finite number of digits (15) used in the calculations (round-off errors). Before 
the saturation, the error is evanescent (exponential decay with the number of collocation points), which is typical of 
spectral methods. 
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B. Kerr black hole 



In this section we consider a single rotating black hole by setting /? / 0. Let us mention that, since the Kerr solution 
is known to diverge from conformal flatness (see e.g. |^^), we will no be able to recover exactly the Kerr metric. In 
other words the obtained solution is expected to violate some of the 5 Einstein equations we decided to ignore. 

So we solve Eqs. (^), (|^) and (||) with boundary conditions (|Tl|) (^2|) and (|l^) on one single sphere. The values 
at infinity are chosen to maintain asymptotical flatness by using Eq. (|^) and (10). The two parameters of our 
rotating black hole are the radius of the throat S and the rotation velocity fl. The total mass M and and angular 
momentum J are computed at the end of the iteration. 

Initialy the values of N and are set to those of a Schwarzschild black hole and the shift is set to zero. Relaxation 
is used for all the fields with a parameter A = 0.5. As for the Schwarzschild computation, we use two shells with the 
same number Nr x Nq x N^p of collocation points in the two shells and in the external compactified domain. The 
iteration is stopped when the relative difference between the shifts obtained at two consecutive steps is smaller than 
5/3 = 10-10. 

Before comparing the obtained solution to the Kerr metric we perform some self-consistency checks, by varying the 
number of coefficients of the spectral expansion. First of all, we need to verify that the regularization function applied 
to the shift by means of Eq. (39) has gone to zero at the end of the computation. Figure H shows that, for various 
values of the Kerr parameter J/M^, the relative norm of the regularization function decreases very fast, as the number 
of coefficients increases. The saturation value of IQ-^^ is due to the criterium we choose to stop the computation 
Sf3 = 10~^'^. Had it been conducted for a greater number of steps, the saturation level of the double precision would 
have been reached. Figure ^ enables us to say that the shift solution of (|4|) fulfills the regularity conditions (^2|) for 
the extrinsic curvature tensor. Let us mention, that the fact that the conformal approximation is not valid, does not 
prevent the correction function /3cor from going to zero. 
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FIG. 3. Relative norm of the regularization function given by Eq. (^^ with respect to the Kerr parameter J/M^, for 
various numbers Nr x Ne x A*'^ of collocation points. 

As seen in Paper I, the total angular momentum can be calculated in two different ways, using a surface integral 
at infinity: 



(49) 



(where m :— d/dip) or an integral on the throat: 



(50) 



where dSi denotes the surface element with respect to the flat metric f. 
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The two results will coincide if and only if the momentum constraint 



D 



(51) 



has been accurately solved in all the space. This is a rather strong test for the obtained value of A^^ . Figure ^ shows 
that the relative difference between the two results rapidly tends to zero, as the number of coefficients increases. The 
same saturation level as in Fig. ^ is observed. 

The last self-consistency check is to verify the virial theorem considered in Sec. ^ In other words we wish to check 
if the ADM and Komar masses are identical, which should be the case for a Kerr black hole. We plotted the relative 
difference between these two masses, for various numbers of collocation points and rotation velocities in Fig. ^. Once 
more this difference rapidly tends to zero as the number of coefficient increases. Contrary to the case of two black 
holes, the angular velocity J7 is not constrained by the virial theorem, reffecting the fact that an isolated black hole 
can rotate at any velocity (smaller than the one of an extreme Kerr black hole). 
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for the relative difference between the angular momentum calculated by means of Eq. (E9[) and 
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FIG. 5. Same as Fig. ^ for the relative error on the virial theorem. 

To end with a single rotating black hole, we check how far the numerical solution is from an exact, analytically 
given, Kerr black hole. Given the ADM mass M and the reduced angular momentum a =^ J /M , an exact Kerr metric 
in quasi-isotropic coordinates would take the form 

ds^ = -NL„dt^ + Bl,y sin2 e [d^ - N^,,,dtf + A^,,, {dr^ + r'dd^) , (52) 

with A'^Korr, -^Korr' ^Kcrr and i?Korr known functions. It obviously differs from asymptotical flatness because A =/= B 
for a 7^ 0. So we define a pseudo-Kerr metric by setting B = A, which gives 



sin^e {dip -N^^^.^dtf + {dr^+r^d9^) , (53) 



where ^'xcrr — ^Kcn- After a numerical calculation, we compute the global parameters M and a, calculate the 
functions iVKerr, -^Korr ^^'^ ^Kcrr and Compare them to the ones that have been calculated numerically. Note that 
N'-'' := — ri. The coefficients of the pseudo-Kerr metric are given by 

" ^ E2 (^2 + ^2) + 2a2SMi?sin2 ^ ^ 

^, 2M iM^ + areas' 9 (AP ~ a^) M (hP - a^f , , 



E(ii'2 + a2) + 2a2Afi?sin2 0' 



where 



A'P - a2 

R:=r+ — + M (57) 

4r 

S := i?2 + cos^ 6*. (58) 

Those analytical functions are then compared with that obtained numerically (see Fig. |^). As expected the 
difference between the fields is not zero and it increases with fi, reflecting the fact that a Kerr black hole deviates 
more and more from conformal flatness as J increases. 

To summarize the results about a single rotating throat, we are confident in the fact that the Eqs. (||), (Q) and ^ 
have been successfully and accurately solved, with the appropriate boundary conditions. On the other hand we do 
not claim to recover the exact Kerr metric, for this latter differs from conformal flatness. 



ie-oi 




FIG. 6. Relative difference between the pseudo-Kerr quantities defined by Eqs. (|54|)-(|56[) and the numerically calculated 
ones with respect to the angular velocity. The computation has been performed with Nr x Ne x A^ = 25 x 17 x 16. The circles 
denote the differences on A, the squares on and the diamonds on N"^ :— P'^ — Q. 
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C. The Misner-Lindquist solution 



Misner ^ and Lindquist |^ have found the conformal factor ^! of two black holes in the static case, i.e. when /3 = 
(see also Ref. [|4| and Appendices A and B of Ref. [Q). In such a case the equation for "if is only 

A^- = 0, (59) 

which was to be solved using boundary conditions (^^ and (|l^). In the case of identical black holes, that is for two 
throats having the same radius a, the solution is analytical and does only depend on the separation parameter 

d , , 

D:=- , 60 
a 

d being the coordinate distance between the centers of the throats. To check if our scheme enables us to recover such 
a solution, we solve Eq. ( |5^ ) with the boundary conditions (|l^) and (|l3|) . We then compute the ADM mass by means 
of the formula (see Paper I) 

M = i &^dSi (61) 
27r Joe 

and compare the result to the analytical value given by a series in Lindquist article j^]. 

Let us mention that, even if Eq. ( p9[ ) is a linear equation (the source is zero), the problem has to be solved by 
iteration because of our method for setting the boundary condition (|l3|). The computation has been conducted with 
a relaxation parameter A — 0.5 and until a convergence of S'if — 10" has been attained. The comparison between 
the analytical and calculated ADM masses is plotted on Fig. ^ for various values of the separation parameter D and 
various numbers of coefficients. The agreement is very good for every value of D. As for the Kerr black hole, when 
the number of coefficients increases, we attain the saturation level of a few 10"^" is due to the threshold chosen for 
stopping the calculation. This test makes us confident about the iterative scheme used to impose boundary conditions 
onto the two throats 5*1 and 5*2. 
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FIG. 7. Relative difference between the calculated and analytical ADM mass for the Misner-Lindquist solution, 
computation has been performed with various number of coefficients A'^,. x Ne x N^. 



The 



To go a a bit further and check the decomposition of the sources into two parts, presented in Sec. II A , we wish 



to consider a test problem with a source different from zero. To do so we consider the Misner-Lindquist problem but 
decide to solve for the logarithm of $ = In^P. The equation for <I> is 



(62) 
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and it must be solved with the following boundary conditions 

$ — > when r — > cxd 



dri 



1 

2^ 



and 



1 

2^' 



s, -"1 dr2 
The source of the equation for $ containing $ itself, it is split as described in Sec. II A by 

a being 1 or 2. At the end of a computation, we compute the ADM mass by using 

1 



M 



2tt 



D'<^dSi 



(63) 
(64) 

(65) 
(66) 



and compare it with the analytical value. The computation used a relaxation parameter A — 0.5 and has been stopped 
when the threshold 5^ = 10~^ has been reached. 
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FIG. 8. Relative difference between the calculated and analytical ADM mass for the Misner-Lindquist solution calculated 
using $ = In^", with respect to the number of coefficients in tp {Ne = A^'^^ + 1 and Nr = 2N^ + 1). The separation parameter 
is D = 10. 

Figure || shows the resulting relative error estimated by means of the ADM mass for = 10 and various numbers of 
coefficients. The convergence is evanescent, i.e. it is exponential as the number of coefficients increases. Unfortunately, 
this convergence is much slower than when the solution was computed using ^ and Eq. (^9|). This comes from the 
very nature of the source of Eq. (|5|). The part of the equation split on the coordinates centered around throat 
1 is the sum of two terms. The first one —Dk^iD'^^i is centered around hole 1 and well described by spherical 
coordinates associated with this hole. We do not expect any problems with this term. The other term is —Di^^iD^^2 
and contains a part that is centered around hole 2. Describing this part using spherical coordinates around hole 1 is 
much more tricky and a great number of coefhcients, especially in (p, is necessary to do it accurately. It is the presence 
of such a component at the location of the other hole that makes the convergence of the calculation much slower in 
this case. Of course, we expect to recover this effect in the calculation of orbiting black holes. 



IV. SEQUENCE OF EQUAL MASS COROTATING BLACK HOLES IN CIRCULAR ORBIT 

A. Numerical procedure 



In this section we concentrate on equal mass black holes. The only parameter is the ratio D between the distance 
of the centers of the holes and the radius of the throats [see Eq. (60)]. We solve Eqs. (||), (Eh and (0), with values at 
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infinity given by (||), ^ and jl^ ) and with boundary conditions on the horizons by (O), ( p^ and (|l3|). We solve for 
various values of fl and choose for solution the only value that fulfills the condition ([141). It turns out that this process 
uniquely determines the angular velocity. Let us call iltrue the only value that equals the ADM and the Komar-like 
masses. It happens that 



if il < ritiuo then M-^ 



< Ma 



DM 



• if > r^truc then AfKomar > MadM • 

The fact that — iltruo has always the same sign than MKomar — -^adm enables us to implement a very efficient 
procedure to determine the orbital velocity. It is found as the zero of the function Mxcmar (^) — -^adm {^) by means 
of a secant method. This is illustrated by Fig. ^, which shows the value (Madm — AfKomar) /^Komar for various 
values of fi, with respect to the step of the iterative procedure. Those calculations have been performed for D = 16. 
The solid line denotes ritruc, the only value of fl for which (Madm — Mxcmar) /-MKomar converges to 0. 
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FIG. 9. Value of (Madm — MKomar) /MKomar with respect to the step of the iteration, for D — 16 and for various values of 
O. The solid line denotes f2truo, the short-dashed line Q = 0.95 Sltruc and the long-dashed line Q, = l.OSiltruc. 

The computations have been done either in low resolution with Nr x Ng x = 21 x 17 x 16 coefficients in each 
domain or in high resolution with Nr x Ng x = 33 x 21 x 20 coefficients in each domain. All the computation 
used a relaxation parameter A = 0.5. We solve first for the static case 51 = and use that solution as initial guess. 
For each value of fl, the computation is stopped for a relative change on the shift vector as small as 6/3 = 10~^ (resp. 
6/3 — 10^^) for the high resolution (resp. low resolution) between two consecutive steps. The secant procedure for the 
determination of the angular velocity has been conducted until |(Madm — MKomar) /MKomar I < 10^^ (resp. < 10^'') 
for the high resolution (resp. low resolution), which gives a precision on Sltruo of the order of 10^^ (resp. 10^"^). 



B. Tests 



1 . Check of the momentum constraint 



As discussed in Sec. II C 1 , we have to slightly modify the shift vector to ensure both the regularity of the extrinsic 
curvature on the throats and the invariance of the shift under the inversion isometry. This modification of the shift, 
via the addition of the correction function ficor, results in a slight violation of the momentum constraint Eq. (^. A 
good way to measure the magnitude of this violation is to check whether the total angular momentum J has the same 
value when calculated by surface integrals at infinity or on the throats. Indeed, as for the Kerr black hole, it has been 
shown in Paper I that J can be given either by one of the two following integrals: 
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J 



1 

8^ 



Si 



Stt 



S2 



(67) 
(68) 



Any difference between those two formulas would reflect the fact that the momentum constraint (|51| ) is not exactly 
fulfilled. 
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FIG. 10. Relative difference between the regularized shift and the exact solution of (|^) (circles) and relative difference 
between J calculated by means of (^^ and (^^ (squares). The filled symbols and the solid line denote the high resolution and 
the empty symbols and the dashed line the low resolution. 

We have plotted the relative difference 6 J/ J between the two integrals ( |67|) and (|6^ ) in Fig. |l^ as a function on 
the separation between the two holes. Also shown on the same figure is the relative norm of the correction function 
|/3cor|/|/3|- The correlation between the two curves shows that the error on the momentum constraint arises from 
the introduction of the correction function on the shift. It is also clear from Fig. |o| that increasing the number of 
coefficients in the spectral method does not make the correction function tend to zero. This means that the error in 
the momentum constraint come rather from the method (necessity to regularize the shift vector) than from some lack 
of numerical precision. 

As discussed in Paper I, we had to regularize the shift vector because Eq. (^Tj) is not enforced in our scheme. It 
has been argued recently by Cook that if one reformulates the problem by assuming that the helical vector £ is 
not an exact Killing vector, but only an approximate one — as it is in reality — then the only freely specifiable part 
of the extrinsic curvature, as initial data, is (^), not (E). This means that the relation ( ^l|) 
curvature and the shift is not as robust as the relation (r^ 

However, we see from Fig. ^that at the innermost sta 
the error are very small: 



between the extrinsic 



)le circular orbit, which is located at D = 17 (cf. Sec. IVC), 



(5J/J = 2 IQ-^ 
|/3cor| = 8 10-^1/31 



(69) 
(70) 



The 6 J/ J error estimator maximizes the error on the momentum constraint because it integrates it in all space. Thus 
we conclude that momentum constraint is satisfied in our numerical results with a precision of the order 1%. 



2. Check of the Smarr formula 



A good check of the global error in the numerical solution is the generalized Smarr formula derived in Paper I : 
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M - 2nj = — - 

An 



4:71 



(71) 



S2 



For any computation, one gets M, Q and can compute the r.h.s. of Eq. ( |71| ) and use that equation to derive the 
value of J that fulfiUs the Smarr formula. That value is then compared to the ones calculated using Eqs. (67) and 
(|68|). The comparison is plotted in Fig. |l^ for 



|). The comparison is plotted in Fig. |ll| for the two different resolutions. It turns out that the angular momentum 
calculated at infinity is better in fulfilling the Smarr formula than the one calculated on the throats by an order of 
magnitude and that the precision is better than 5 x 10"'^. So, for all following purposes, we will use the value of J 
given by Eq. (|67|). 
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FIG. 11. Relative error on the generalized Smarr formula (|7l|). The circles denote the error obtained using J calculated at 
infinity [Eq. (|67|)] and the squares that obtained when evaluating J on the throats [Eq. (|6^]. The filled symbols and the solid 
line denote the high resolution and the empty symbols and the dashed line the low resolution. 



3. Check of Kepler law at large separation 

The next thing one wishes to test is the value of fi, obtained from the virial criterium (p^). In Newtonian gravity, 
two points particles on circular orbits obey the following relation, which is equivalent to Kepler's third law: 

where M is the total mass, J the total angular momentum and the orbital velocity. For large separations of the 
two throats we expect to recover this relation. Therefore, for every value of D, we evaluate 



4jr!i/3 

1 

and check if / tends to 1 when D ^ 00. 
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FIG. 12. Value of / = 4J {Q/M^Y^'^ (1 ow resolution runs) with respect to the separation parameter D. The horizontal 
dashed line corresponds to the value predicted by Kepler's third law. 

The value of / is plotted in Fig. |l^ with respect to the distance parameter D. As expected, for large values of D, 
it tends to 1, implying that for large separations the system behaves like two point particles in Keplerian motion. 



C. Evolutionary sequence 



Let us first present some figures about the metric fields. Figure 13 shows the total lapse function N, conformal 
factor 5* and the shift vector f3 and Fig. |lj the components A^^, A^^ and A^^ of the extrinsic curvature tensor. 
All those plots are cross-section in the orbital plane z = and the coordinate system is a Cartesian one centered at 
the middle of the centers of the throats. The computation has been done using the high resolution. The separation 
parameter is D = 17. As it will be seen later, this separation corresponds to the turning point in the energy and 
angular momentum curves. 



Lapse function (Z=0) 



Conformal factor (Z=0) 



Slnft vector (Z=0) 




FIG. 13. Isocontour of the lapse function A'^ and of the conformal factor <!/ and plot of the shift vector /3, for D = 17, in the 
orbital plane z — 0. The computation has been done using the high resolution. The thick solid lines denote the surfaces of the 
throats. The kilometer scale of the axis corresponds to an ADM mass of 31.8M0. 
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FIG. 14. Isocontour of the extrinsic curvature tensor for D = 17 in the orbital plane z = 0. The solid (dashed) lines denote 
positive (negative) values. The thick solid lines denote the surfaces of the throats. The computation has been done using the 
high resolution. The kilometer scale of the axis corresponds to an ADM mass of 31.8Mq. 

In the previous section, the only parameter we considered was the dimensionless separation parameter D. But 
there also exists a scaling factor. Suppose that all the distances in the computation are multiplied by some factor a. 
Another solution with the same value of D will be obtained, the global quantities being rescaled as 

M„ = aMi (74) 
J„ = aVi (75) 

= — , (76) 
a 

where Mi, Ji and Qi are the values before rescaling and Ma, Ja and Cla the values after the rescaling. 

Consider a physical configuration corresponding to a value D (n) of the separation parameter, with global quan- 
tities M (n), J {n) and Q (n). This system will evolve due to the emission of gravitational radiation. A subsequent 
configuration n + 1 will have _D (n -j- 1) < D (n). But what scaling factor a should be applied to the configuration 
calculated for D (n + 1) to ensure that it represents the same physical system as before ? In other word, a physical 
sequence is a one parameter (the separation) family of configurations and we have to impose another condition to 
determine the scaling factor associated with each value oi D. In the case of binary neutron stars the condition is 
obtained by imposing that the number of baryons is conserved (see e.g. Ref. Q). This cannot be extended to the 
black holes case since no matter is present. We chose instead to define a sequence by requiring that the loss of energy 
(ADM mass) dAI and angular momentum dJ due to gravitational wave emission are related by 



dM 
U 



= n. (77) 

sequence 



This relation is exact at least when one considers only the quadrupole formula (see e.g. page 478 of Ref. It 



turns out that it is also well verified for sequences of binary neutron stars [g7 28|. So Eq. (77) should hold rather well 
for corotating black holes. 

The scaling factor a associated with the separation parameter D {n+ 1) can be computed from the global values 
at separation D (n) and the unsealed values at separation D {n+ 1) as the solution the third degree equation 

M(n)-aMi(n + l) 1 f^, , , n,{n+l) \ 

J{n)-aU,{n+l) = H^"^ ^ J ' ^ ^ 

which is a first-order translation of Eq. (|7^). To present the results, we define the following dimensionless quantities 

M 

M:.- (79) 
J := ±, (80) 



19 



n := Man (81) 
1:=—, (82) 

where I is the proper separation of the holes, defined as the geometrical distance between the throats along the axis 
joining their centers. Mg is some arbitrary mass used for normalization purpose. It is often convenient to choose Mg 
to be the total mass of the system when the two holes are infinitely separated, i.e. the ADM mass when D ^ oo. 
Unlike other methods, this value is not an input parameter of our calculation. It can only be obtained by constructing 
a sequence until very large values of D, which would impose to calculate a great number of configurations. However, 
as will been seen further, the system will exhibit turning point in the total energy and angular momentum, thereafter 
assumed to be the signature of an innermost stable circular orbit (thereafter ISCO). We chose Mg to be the total 
ADM mass of the system at that point : 

Mo := MadmIiscO' (^^) 

so that AI is 1 at the location of the ISCO. 

The values of the dimensionless quantities f2, J, M and I along the sequence are given by Table |, for the high 
resolution. 

TABLE I. Values of dimensionless quantities along a sequence of corotating black holes obtained using the high resolution. 
The bold line denotes the values at the location of the ISCO. 
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FIG. 15. M with respect to J along a sequence. The filled symbols and solid line denote the high resolution and the empty 
symbols and the dashed line the low one. 
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FIG. 16. M with respect to Q. along a sequence. The filled symbols and solid line denote the high resolution and the empty 
symbols and the dashed line the low one. 
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FIG. 17. J with respect to Q. along a sequence. The filled symbols and solid line denote the high resolution and the empty 
symbols and the dashed line the low one. 
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FIG. 18. I with respect to S7 along a sequence, 
symbols and the dashed line the low one. 



The filled symbols and solid hne denote the high resolution and the empty 



Figures |l^, ^ and |l^ show the values of the dimensionless quantities along a sequence. The calculation has been 
performed with the high and low resolutions and for values of the parameter D ranging from 40 to 11. As previously 
mentioned, the sequence exhibits a minimum of J and M as the throats become closer, thereafter interpreted as the 
signature of an innermost stable circular orbit (ISCO) But at this point, we have to be cautious. Indeed, the 

relative variation of M and J a long a sequence is rather small, and comparable to the precision estimated by means 



of the Smarr formula (see Sec. IV B ). The exact location of the ISCO being very dependent on those small effects, 
we do not claim to have very precisely determined it. The following results should be confirmed with more precise 
calculations. 
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Another important quantity is the area of the black hole horizons which relates to the irreducible mass (see 
also Box 33.4 of |^^). As discussed in Sec. II. B. 6 of Paper I, in our case the apparent horizons coincide with the two 
throats. We therefore define the dimensionless irreducible mass by 




where Aa (a = 1, 2) denotes the area of the throat a, calculated according to the formula 



Aa = 



(84) 



(85) 
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FIG. 19. Relative change of Mir along the sequence, with respect to the orbital velocity $7. The filled symbols and the solid 
line denote the high resolution and the empty symbols and the dashed line the low resolution. 



Figure shows the relative change of M-^ along the sequence. It exhibits a slight increase, but its variation is very 
small. It appears that, along the overall sequence, the variation is smaller than 10~^. The precision of our code being 
of that order, this result is compatible with the fact that M-„ is constant. In other word it shows that imposing the 
condition ( |77| ) is equivalent to imposing that the irreducible mass is constant along a sequence. This constitutes in 
fact a very good test of our procedure. Indeed Friedman, Uryu & Shibata have recently established the first law 
of binary black hole thermodynamics: 



dM = ridJ + KidAi + K2dA2 , 



(86) 



where ni and K2 are two constants, representing the black holes surface gravity. For identical black holes (ki = K2 
and dAi = dA2), the above relation implies 



dM = ndj 



dAa =0 (a = 1, 2) . 



(87) 



Hence the area of each black hole must be conserved during the evolution. In future works, this last criterium could 
be used to define a sequence, instead of the relation ([ttI). 

We choose an average value of the irreducible mass Mi^ = 1.0173 and we define then the binding energy of the 
system at the ISCO by 

^bliscO = 1 -^ir. 



the dimensionless total energy being equal to 1 at the location of the turning point. 
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The values of the dimensionless quantities fl, J, E\y and I at the ISCO are given in Table || and compared with 
the results from other approaches (see for a review). 3-PN EOB stands for the third order post-Newtonian 
Effective One Body method for non-spinning black holes , with two values of the 3-PN parameter : ujs — Q and 
LOs = —9.34. 3-PN j-method stands for third order post-Newtonian j-method of |^^. Puncture denotes the results 
from the puncture method in the case of non-spinning black holes |34 and Conformal imag. the conformal imaging 
approach with various values of the individual spins for rotating black holes (111 ■ definition M = 1 at the location 
of the ISCO for all the methods. The results from the different methods are also plotted in Fig. pO|. 
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FIG. 20. Values of i?b and J with respect to Q, at the ISCO for different methods, including ours with high and low resolution. 
The references to previous studies are as follows: Damour et al. 2000: |Q, Pfeiffer et al. 2000: and Baumgarte 2000: |Q. 
S denotes the (fixed) spin of the black holes used in various methods. 



TABLE II. Values of dimensionless quantities at the location of the ISCO. Comparison with other works. 



Method 

3-PN EOB a;, = S = H 
3-PN EOB ujs = -9.34 S^Q\ 
3-PN j-method c^^ = -9.34 S = 
Puncture S = |@ 
Conformal imag. 5 = p4| 
Conformal imag. S = 0.08 
Conformal imag. S = 0.17 in! 
This work (high res.) 
This work (low res.) 
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J 




I 


0.0868 


0.847 


-0.0170 


not given 


0.0722 


0.877 


-0.0152 


not given 


0.0731 


0.877 


-0.0153 


not given 


0.176 


0.773 


-0.0235 


4.913 


0.162 


0.779 


-0.0230 


5.054 


0.182 


0.799 


-0.0250 


4.705 


0.229 


0.820 


-0.0279 


4.040 


0.101 


0.869 


-0.0173 


6.606 


0.105 


0.867 


-0.0173 


6.450 
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Figure shows explicitly that the present results are in much better agreement with post-Newtonian calculations 
than with other numerical works. Note that the post-Newtonian point plotted on Fig. ^ corresponds to the value 
of the (previously ambiguous) 3-PN "static" parameter uig- This is indeed the value recently determined by Damour 
et al. pSfl . 

But let us point out that it is rather difficult to compare precisely our results with the other works. The main 
problem comes from the fact that all those methods use individual spins of the black holes as input parameters. In the 
present paper we impose corotation, that is that the throats are spinning at the orbital velocity. The only value that 
can be computed is the total angular momentum J and, in general relativity, it cannot be split into orbital and spins 
parts in a invariant way. However, from the results of Pfeiffer et al. jl^] one can see that increasing the spins of the 
black holes make the values f2, J and —E^ at the ISCO greater. Taking rotation into account in the post-Newtonian 
methods [p6| will probably make the orbital velocity and the binding energy at the ISCO match even better with our 



values. Work is under progress to compare with corotating post-Newtonian results |37 



So, it appears that our results match pretty well with post-Newtonian methods. This is the most striking conclusion 
from our study. The difference between numerical and post-Newtonian results have often been imputed mostly to 
the conformal flatness approximation (see ||29tl ). The fact that our result, using conformal flatness, is in much better 
agreement with FN calculations than other numerical works, makes us believe that the main worry of both conformal 
imaging and puncture methods lies elsewhere, possibly in the determination of Q. Indeed, it is very unlikely that 
the orbits and so orbital velocity can be properly computed by solving only for the four constraint equations. Time 
should be involved at some level and one should take other Einstein equations into account, as we have done here. 



V. CONCLUSIONS 



The present work should be seen as a first step in trying to give some new insight to the binary black holes problem. 
The basic idea is to extend the numerical treatment beyond the resolution of the four constraint equations within a 
3-dimensional spacelike surface. This is achieved by reintroducing time in the problem to deal with a 4-dimensional 
spacetime. The orbits are well defined by imposing the existence of a helical Killing vector and the orbital velocity 
is found as the only value that equals the ADM and the Komar-like masses, a requirement which is equivalent to 
the virial theorem. According to us those are the two most important features of our method. The approximation 
of conformal flatness for the 3-metric has only been used for simplicity. Sooner or later this problem will have to 
be solved using a general spatial metric and outgoing waves boundary conditions at large distances. The use of the 
inversion isometry to derive boundary conditions on the throats is also a weak assumption. In the future, it would be 
interesting to change the boundary conditions on the fields in order to investigate their influence on the results (see 
e.g. for an alternative choice). Besides, changing the boundary conditions on the shift vector should enable us to 
describe other states of rotation of the black holes, as has been recently proposed by Cook p^ . 

The numerical schemes are basically the same as those which have been previously successfully applied to binary 
neutron stars configurations [Q. They have been extended to solve elliptic equations with non-trivial boundary 
conditions imposed on two throats and exact boundary conditions at infinity. Those techniques passed numerous 
tests and recover the Schwarzschild and Kerr solutions as well as the Misner-Lindquist one for two static black holes 
. A technical problem lies in the great number of coefficients needed to accurately describe the part of the sources 
located around the companion hole. This effect causes some lack of precision. But we can estimate the error it 
generates by varying the number of coefficients, and comparing the results. This is what we have done here, using 
21 X 17 X 16 coefficients in each of the 12 domains for the low resolution computations and 33 x 21 x 20 coefficients 
for the high resolution ones. The accuracy, estimated from the generalized Smarr formula, is below 1%. 

Another issue is the slight violation of the momentum constraint which arises from the necessity to regularize the 
shift vector. We have found that the modification of the shift vector with respect to the vector which satisfies the 
momentum constraint (jj) is below 10~^, and that the error it induces in the momentum constraint equation is of 
the order 1%. In view of the other approximations performed in this work, especially the conformal flatness of the 
3-metric, we find this to be very satisfactory. 

In this article, we have defined a sequence of binary black holes by requiring that the ADM mass decrease is related 
to the angular momentum decrease via dAI = Q dJ. This relation is true for the loss due to gravitational radiation, 
at least when one considers only the quadrupole formula. We have then found that the area of the apparent horizons 
(irreducible mass) is constant along the sequence, in agreement with the first law of binary black holes thermodynamics 
recently derived by Friedman et al. fs^ . 

The location of the ISCO has been obtained and compared with the results from other methods p3| , p^ , pT[ . It 
turns out that our results match the 3-PN methods much better than previous numerical works. The differences 
between numerical studies and 3-PN approximations have often been explained by the use of the conformal flatness 
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approximation in the numerical calculations [ p3[ . It seems to us that this is not the main explanation, for we are 
using this approximation. It certainly arises instead from the way is determined. 

Another natural extension of this work is to use the obtained configurations as initial data for binary black holes 
evolution codes (see for a review and Rcfs. [^-|4ll for recent results). Initial data files containing the result 
of the present work are publically available on the CVS repository of the European Union Network on Sources of 
Gravitational Radiation [ff2| . Extraction of the wave- forms from a sequence would also be an interesting application 
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